/*
 *  sparsematrixslu.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos García Orden on 01/06/10.
 *
 */

#include "sparsematrixslu.h"

sparseMatrixSLU::sparseMatrixSLU ()
:	sparseMatrix("slu"),
	A(),
	L(),
	U(),
	B(),
	X(),
	lwork(0),
	work(),
	options(),
	stat(),
	rhsb(),
	rhsx(),
	perm_r(),
	perm_c(),
	etree(),
	equed(),
	R(),
	C(),
	rpg(),
	rcond(),
	ferr(),
	berr(),
	mem_usage(),
	info()
{
	setTypeName("slu");
}

// ***

sparseMatrixSLU::sparseMatrixSLU (int nrow, int ncol)
:	sparseMatrix("slu",nrow, ncol),
	A(),
	L(),
	U(),
	B(),
	X(),
	lwork(0),
	work(),
	options(),
	stat(),
	rhsb(),
	rhsx(),
	perm_r(),
	perm_c(),
	etree(),
	equed(),
	R(),
	C(),
	rpg(),
	rcond(),
	ferr(),
	berr(),
	mem_usage(),
	info()
{
	setTypeName("slu");
}

// ***

sparseMatrixSLU::~sparseMatrixSLU ()
{
	if(&stat!=0) StatFree(&stat);
	
    SUPERLU_FREE (etree);
	SUPERLU_FREE (perm_r);
	SUPERLU_FREE (perm_c);
	SUPERLU_FREE (R);
	SUPERLU_FREE (C);
	SUPERLU_FREE (ferr);
	SUPERLU_FREE (berr);
	Destroy_CompCol_Matrix(&A);
	
    Destroy_Dense_Matrix(&B);
	Destroy_Dense_Matrix(&X);
   	if ( lwork >= 0 )
  	{		/* Deallocate storage associated with L and U. */
		Destroy_SuperNode_Matrix(&L);
		Destroy_CompCol_Matrix(&U);
	}	
}

// ***

void sparseMatrixSLU::initMatrices(void)
{
	int i;
	
	if(nnz<2000) lwork = 0;
	else lwork = nnz*450L;

//cout << "Maximum value for int: " << numeric_limits<long>::max() << endl;	
	cout << "++++++++++++++++++ NNZ = " << nnz << endl;
	cout << "++++++++++++++++++ LWORK = " << lwork/1e9 << " GB " << endl;
	cout << "++++++++++++++++++ LWORK = " << lwork << "  " << endl;
	
	set_default_options(&options);
	options.ColPerm = MMD_AT_PLUS_A;
	
	if ( lwork > 0 ) {
		work = SUPERLU_MALLOC(lwork);
		if ( !work ) {
			ABORT("DLINSOLX: cannot allocate work[]");
		}
    }
	
    dCreate_CompCol_Matrix(&A, n, n, nnz, values, ia, ja, SLU_NC, SLU_D, SLU_GE);

	if ( !(perm_c = intMalloc(n)) ) ABORT("Malloc fails for perm_c[].");
	if ( !(perm_r = intMalloc(n)) ) ABORT("Malloc fails for perm_r[].");
	if ( !(etree = intMalloc(n)) ) ABORT("Malloc fails for etree[].");
	if ( !(R = (double *) SUPERLU_MALLOC(A.nrow * sizeof(double))) ) 
		ABORT("SUPERLU_MALLOC fails for R[].");
	if ( !(C = (double *) SUPERLU_MALLOC(A.ncol * sizeof(double))) )
		ABORT("SUPERLU_MALLOC fails for C[].");

	if ( !(rhsx = doubleMalloc(n * nrhs)) ) ABORT("Malloc fails for rhsx[].");
	if ( !(rhsb = doubleMalloc(n * nrhs )) ) ABORT("Malloc fails for rhsb[].");
	
	if ( !(ferr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )	
		ABORT("SUPERLU_MALLOC fails for ferr[].");
	
	if ( !(berr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )
		ABORT("SUPERLU_MALLOC fails for berr[].");
	
	dCreate_Dense_Matrix(&B, n, nrhs, rhsb, n, SLU_DN, SLU_D, SLU_GE);
	dCreate_Dense_Matrix(&X, n, nrhs, rhsx, n, SLU_DN, SLU_D, SLU_GE);
	
	//cout << "n = " << n  << endl;
	//initVectors(0.005*n);
}

// ***

void sparseMatrixSLU::factorize(int facMode)
{
	((NCformat*) A.Store)->nzval = values;


	if(facMode==0)	options.Fact = DOFACT;	// factorize from scratch
	else	options.Fact = SamePattern;

//	options.Fact = DOFACT;

	StatInit(&stat);
	
	/* ONLY PERFORM THE LU DECOMPOSITION */
    B.ncol = 0;  /* Indicate not to solve the system */

	//dPrint_CompCol_Matrix("A", &A);
	
	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B, &X, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
//	dPrint_CompCol_Matrix("U", &U); 
//	dPrint_SuperNode_Matrix("L", &L);

//	printf("LU factorization: dgssvx() returns info %d\n", info);
	
/*
	
	if ( info == 0 || info == n+1 )
	{
		if ( options.PivotGrowth ) printf("Recip. pivot growth = %e\n", rpg);
		if ( options.ConditionNumber )
			printf("Recip. condition number = %e\n", rcond);
		
		SCformat *Lstore;
		NCformat *Ustore;
		Lstore = (SCformat *) L.Store;
		Ustore = (NCformat *) U.Store;
		printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
		printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
		printf("No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz - n);
		printf("FILL ratio = %.1f\n", (float)(Lstore->nnz + Ustore->nnz - n)/nnz);
		
		printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
			   mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
		fflush(stdout);
	} else if ( info > 0 && lwork == -1 ) {
		printf("** Estimated memory: %d bytes\n", info - n);
	}

	if ( options.PrintStat ) StatPrint(&stat);
*/
	StatFree(&stat);
	
	options.Fact = FACTORED;
}

// ***

void sparseMatrixSLU::subsSolve(vector<double>& b_in)
{
	int i;
	
	// Fill right hand side vector
	for (i = 0; i < b_in.size(); ++i)	rhsb[i] = b_in[i];
	
	((DNformat*) B.Store)->nzval = rhsb;
	
	// Solve, provided matrix is already factorized
	options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
	B.ncol = nrhs;  /* Set the number of right-hand side */
	
	StatInit(&stat);
	
	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B, &X, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
	StatFree(&stat);
	
	// Get solution
	rhsx = (double*) ((DNformat*) X.Store)->nzval;
    for (i = 0; i < n; ++i) b_in[i] = rhsx[i];
}

// ***
//Arma
/*
void sparseMatrixSLU::subsSolve(rowvec & b_in)
{
}
*/
// ***

void sparseMatrixSLU::subsSolve(RowVectorXd & b_in)
{
	int i;
	
	// Fill right hand side vector
	for (i = 0; i < b_in.size(); ++i)	rhsb[i] = b_in(i);
	
	((DNformat*) B.Store)->nzval = rhsb;
	
	// Solve, provided matrix is already factorized
	options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
	B.ncol = nrhs;  /* Set the number of right-hand side */
	
	StatInit(&stat);
	
	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B, &X, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
	StatFree(&stat);
	
	// Get solution
	rhsx = (double*) ((DNformat*) X.Store)->nzval;
    for (i = 0; i < n; ++i) b_in(i) = rhsx[i];
		
}

// ***

void sparseMatrixSLU::subsSolve(VectorXd & b_in)
{
	int i;
	
	// Fill right hand side vector
	for (i = 0; i < b_in.size(); ++i)	rhsb[i] = b_in(i);
	
	((DNformat*) B.Store)->nzval = rhsb;
	
	// Solve, provided matrix is already factorized
	options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
	B.ncol = nrhs;  /* Set the number of right-hand side */
	
	StatInit(&stat);
	
	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B, &X, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
	StatFree(&stat);
	
	// Get solution
	rhsx = (double*) ((DNformat*) X.Store)->nzval;
    for (i = 0; i < n; ++i) b_in(i) = rhsx[i];
	
}

// ***
/*
void sparseMatrixSLU::subsSolve(vector<vector<double> >& B_in)
{
	int i;
	int j;
	
	int ncolrhs = B_in[0].size();
	
	// Give proper size to B
	initVectors(ncolrhs);
	
	// Fill right hand side vector (from a matrix)
	for(j=0; j<ncolrhs; j++)
	{
		for(i=0; i<B_in.size(); i++)
		{
			rhsb[i+B_in.size()*j] = B_in[i][j];
		}
	}
	
	((DNformat*) B.Store)->nzval = rhsb;
	
	// Solve, provided matrix is already factorized
	options.Fact = FACTORED; // Indicate the factored form of A is supplied.
	B.ncol = nrhs;  //Set the number of right-hand side
	
	StatInit(&stat);
	
	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B, &X, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
	StatFree(&stat);
	
	// Get solution
	rhsx = (double*) ((DNformat*) X.Store)->nzval;
	for(i=0; i<B_in.size(); i++)
	{
		for(j=0; j<ncolrhs; j++)
		{
			B_in[i][j] = rhsx[j + ncolrhs*i];
		}
	}		
}
*/
// ***

void sparseMatrixSLU::subsSolveTrans(vector<vector<double> >& B_in)
{
	int i;
	int j;
	
	int ncolrhs = B_in.size();
	int nrowrhs = B_in[0].size();
	
//	double t1 = MPI::Wtime();
	// Give proper size to B
	initVectors(ncolrhs);
//	double t2 = MPI::Wtime();
//	cout << "En subsSolveTrans el tiempo de reservar memoria es: " << t2-t1 << endl;
	
//	t1 = MPI::Wtime();
	// Fill right hand side vector (from a matrix)
	for(i=0; i<ncolrhs; i++)
	{
		for(j=0; j<nrowrhs; j++)
		{
			rhsb[j+nrowrhs*i] = B_in[i][j];
		}
	}
	
	((DNformat*) B.Store)->nzval = rhsb;

	// Solve, provided matrix is already factorized
	options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
	B.ncol = nrhs;  /* Set the number of right-hand side */
	
//	t2 = MPI::Wtime();
//	cout << "En subsSolveTrans el tiempo de llenar vector rhsb es: " << t2-t1 << endl;	

	//t1 = MPI::Wtime();
	
	StatInit(&stat);

	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B, &X, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
	StatFree(&stat);

//	t2 = MPI::Wtime();
//	cout << "En subsSolveTrans el tiempo de resolver es : " << t2-t1 << endl;	
	
//	t1 = MPI::Wtime();
	// Get solution
	rhsx = (double*) ((DNformat*) X.Store)->nzval;
	for(i=0; i<ncolrhs; i++)
	{
		for(j=0; j<nrowrhs; j++)
		{
			B_in[i][j] = rhsx[j+nrowrhs*i];
		}
	}		
//	t2 = MPI::Wtime();
//	cout << "En subsSolveTrans el tiempo de volver a llenar B_in es : " << t2-t1 << endl;		
}

// ***
/*
void sparseMatrixSLU::subsSolveTrans(double* rb)
{
	cout << "HHHHHooolalllaaa " << endl;
	// Solve, provided matrix is already factorized
	options.Fact = FACTORED;
	//B_in.ncol = nrhs; 
	
	//X = B_in;
	nrhs = 30;

	if ( !(rhsx = doubleMalloc(n * nrhs)) ) ABORT("Malloc fails for rhsx[].");
	
	dCreate_Dense_Matrix(&X, n, nrhs, rhsx, n, SLU_DN, SLU_D, SLU_GE);
	
	StatInit(&stat);
	
	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B_in, &X, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
	StatFree(&stat);
}
*/
// ***

	//Arma
/*
void sparseMatrixSLU::subsSolveTrans(mat & B_in)
{
	int i;
	int j;
	
	int ncolrhs = B_in.n_rows;
	
	// Give proper size to B
	initVectors(ncolrhs);
	
	// Fill right hand side vector (from a matrix)
	for(i=0; i<ncolrhs; i++)
	{
		for(j=0; j<B_in.n_cols; j++)
		{
			rhsb[j+B_in.n_cols*i] = B_in(i,j);
		}
	}
	
	((DNformat*) B.Store)->nzval = rhsb;
	
	// Solve, provided matrix is already factorized
	options.Fact = FACTORED; 
	B.ncol = nrhs;
	
	StatInit(&stat);
	
	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B, &X, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
	StatFree(&stat);
	
	// Get solution
	rhsx = (double*) ((DNformat*) X.Store)->nzval;
	for(i=0; i<ncolrhs; i++)
	{
		for(j=0; j<B_in.n_cols; j++)
		{
			B_in(i,j) = rhsx[j+B_in.n_cols*i];
		}
	}			
}
*/
// ***

void sparseMatrixSLU::subsSolveTrans(MatrixXd & B_in)
{
	int i;
	int j;
	
	int ncolrhs = B_in.rows();
	
	// Give proper size to B
	initVectors(ncolrhs);
	
	// Fill right hand side vector (from a matrix)
	for(i=0; i<ncolrhs; i++)
	{
		for(j=0; j<B_in.cols(); j++)
		{
			rhsb[j+B_in.cols()*i] = B_in(i,j);
		}
	}
	
	((DNformat*) B.Store)->nzval = rhsb;
	
	// Solve, provided matrix is already factorized
	options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
	B.ncol = nrhs;  /* Set the number of right-hand side */
	
	StatInit(&stat);
	
	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B, &X, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
	StatFree(&stat);
	
	// Get solution
	rhsx = (double*) ((DNformat*) X.Store)->nzval;
	for(i=0; i<ncolrhs; i++)
	{
		for(j=0; j<B_in.cols(); j++)
		{
			B_in(i,j) = rhsx[j+B_in.cols()*i];
		}
	}				
}

// ***

void sparseMatrixSLU::subsSolveTrans(SuperMatrix& B_in, SuperMatrix& X_out)
{
//	SuperMatrix Xtemp = Bin;
	
	//options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
	
	//StatInit(&stat);
	
	dgssvx(&options, &A, perm_c, perm_r, etree, equed, R, C,
           &L, &U, work, lwork, &B_in, &X_out, &rpg, &rcond, ferr, berr,
           &mem_usage, &stat, &info);
	
	//StatFree(&stat);	
}

// ***

void sparseMatrixSLU::initVectors(int nvec)
{
	nrhs = nvec;
	if ( !(rhsx = doubleMalloc(n * nrhs)) ) ABORT("Malloc fails for rhsx[].");
	if ( !(rhsb = doubleMalloc(n * nrhs )) ) ABORT("Malloc fails for rhsb[].");
	
	if ( !(ferr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )	
		ABORT("SUPERLU_MALLOC fails for ferr[].");
	
	if ( !(berr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )
		ABORT("SUPERLU_MALLOC fails for berr[].");
					
	dCreate_Dense_Matrix(&B, n, nrhs, rhsb, n, SLU_DN, SLU_D, SLU_GE);
	dCreate_Dense_Matrix(&X, n, nrhs, rhsx, n, SLU_DN, SLU_D, SLU_GE);
}

// ***

void sparseMatrixSLU::copyMatrix(sparseMatrix* b)
{
	sparseMatrixSLU *aux = dynamic_cast<sparseMatrixSLU*> (b);
		
	aux->n = n;
	aux->ia = ia;
	aux->ja = ja;
	aux->nnz = nnz;
	aux->nrhs = nrhs;
	aux->typeName = typeName;
	int i;
	for (i=0; i<nnz; ++i)
	{
//		cout << "value[" << i << "] = " << values[i] << endl;
		aux->values[i] = values[i];
	}
	
	aux->A = A;
	aux->L = L;
	aux->U = U;
	aux->B = B;
	aux->X = X;
	aux->lwork = lwork;
	aux->work = work;
	aux->options = options;
	aux->stat = stat;
	aux->rhsb = rhsb;
	aux->rhsx = rhsx;
	aux->perm_r = perm_r;
	aux->perm_c = perm_c;
	aux->etree = etree;
	aux->equed[0] = equed[0];
	aux->R = R;
	aux->C = C;
	aux->rpg = rpg;
	aux->rcond = rcond;
	aux->ferr = ferr;
	aux->berr = berr;
	aux->mem_usage = mem_usage;
	aux->info = info;
	
	b = dynamic_cast<sparseMatrix*> (aux);
}

// ***

//void sparseMatrixSLU::freeN()
//{}
